Bi-Affinity Filter: A Bilateral Type Filter for Color Images

ABSTRACT

An edge preserving filter that works on the principle of matting affinity allows a better representation of the range filter term in bilateral class filters. The definition of the affinity term can be relaxed to suit different applications. An approximate bi-affinity filter whose output is shown to be very similar to the traditional bilateral filter is defined. The present technique has the added advantage that no color space changes are required and hence an input image can be handled in its original color space. This is a big benefit over the traditional bilateral filter, which needs conversion to perception based spaces, such as CIELAB, to generate results close to the present invention. The full bi-affinity filter preserves very minute details of the input image, and thus permits an enhanced zooming functionality.

BACKGROUND

1. Field of Invention

The present invention pertains to the field of image processing. Morespecifically, it pertains to the field of image smoothing by bilateralfiltering.

2. Description of Related Art

The bilateral filter was originally proposed by Tomasi and Manduchi inBilateral Filtering for Gray and Color Images, ICCV: Proceedings of theSixth International Conference on Computer Vision, page 839, 1998, IEEEComputer Society, herein incorporated in its entirety by reference.Basically, the bilateral filter smooths an image while preserving edgesby means of a nonlinear combination of nearby image values. Theprinciple idea behind such a filtering operation is to combineinformation from the spatial domain as well as from a feature domain. Itcan be represented as

$\begin{matrix}{{h(x)} = {\frac{1}{k(x)}{\sum\limits_{y \in \Omega_{x}}{{f_{s}\left( {x,y} \right)}{g_{r}\left( {{I(x)},{I(y)}} \right)}{I(y)}}}}} & (1)\end{matrix}$

where I and h are the input and output images respectively, x and y arepixel locations over an image grid, Ω_(x) is the neighborhood inducedaround the central pixel x, f_(s)(x,y) measures the spatial affinitybetween pixels at x and y (i.e. the spatial domain, and may be thoughtof as a spatial filter) and g_(r)(I(x),I(y)) denotes thefeature/measurement/photometric affinity (i.e. the feature domain, andmay be thought of as a range filter). Parameter k(x) is thenormalization term given by

$\begin{matrix}{{k(x)} = {\sum\limits_{y \in \Omega_{x}}{{f_{s}\left( {x,y} \right)}{g_{r}\left( {{I(x)},{I(y)}} \right)}}}} & (2)\end{matrix}$

The spatial and range filters (f and g, respectively), are commonly setto be Gaussian filters, as follows:

$\begin{matrix}{{f_{s}\left( {x,y} \right)} = {\exp\left( \frac{- {{x - y}}_{2}^{2}}{2\; \sigma_{s}^{2}} \right)}} & (3) \\{{g_{r}\left( {u,v} \right)} = {\exp\left( \frac{- {{u - v}}_{2}^{2}}{2\; \sigma_{r}^{2}} \right)}} & (4)\end{matrix}$

parameterized by the variances σ_(s) and σ_(r). The range filter, g,penalizes distance in the feature space and hence the filter has aninherent edge preserving property. Due to this property, the bilateralfilter, as an edge-preserving filter, has been one of the most widelyused filtering techniques within computer vision community.

The bilateral filter is a non-linear filter (making it a computationallyintensive operation) and, as such, many researchers have proposedtechniques to decompose the non-linear, bilateral filter into a sum ofseparate one-dimensional filters or similar cascaded representations.Singular value decomposition of the 2D kernel is one such technique.Another proposed technique is an approximation of the bilateral filterby filtering sub-sampled copies of an image with discrete intensitykernels, and recombining the results using linear interpolation.

Recently the run-time of the bilateral filter has been identified as acritical bottleneck, and a few techniques have been proposed that renderthe filtering operation at an almost constant time, albeit with largerspace requirements and behavioral approximations.

The research into improving the performance of the bilateral filterheavily relies on the form of the filter, which is applied in a rangedomain as well as the spatial domain. One method can be entirely brokendown to an approximation of a product of a box filter for smoothing anda polynomial or fourth order Taylor series approximation of a Gaussiankernel. However, when the form of the filter is changed, such methodscannot be applied for the promised increased in speed.

It is an object of the present invention to provide a filter thatprovides similar image smoothing and edge preserving properties as abilateral filter, but which is less constrained by nonlinear combinationof image values.

It is a further object to provide of the present invention to provide afilter that provides similar image smoothing and edge preservingproperties as a bilateral filter, but which lends itself to increasedimplementation optimizations.

It is a further object of the present invention to provide a filter thatnot only provides similar image smoothing and edge preserving propertiesas a bilateral filter, but also provides additional image optimizationproperties.

SUMMARY OF INVENTION

The above objects are met in method of smoothing and enhancing an inputimage by using color line technqiues to determine the color affinity ofpixels within their native color space.

More specifically, the above objects are met in a method of processing adigital input image, comprising: providing the digital input image, I,in physical memory, the input image having color information in a nativecolor space for each image pixel; providing a processing device coupledto the physical memory to access the input image I and create an outputimage, h, by implementing the following steps:

(a) applying the following sequence to the input image I,

${h_{\sigma,\mspace{14mu} ɛ}(x)} = \frac{\sum\limits_{y \in \Omega_{x}}\left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y \in \Omega_{x}}\left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}$

where x and y are pixel locations over an image grid, and Ω_(x) (orequivalently window w) is the neighborhood induced around the centralpixel x; f_(s)(x,y) measures the spatial affinity between pixels at xand y; parameter σ controls an amount of spatial blurring; L_(xy) ^(ε)is a laplacian matrix that provides positive color affinity value fortwo examined pixels, x and y, that have the same color and provides azero color affinity value for pixels with different color; parameter εis a regularization term whose relative weight determines the smoothnessof output image h.

In this method, step (a) is applied to the input image I in its nativecolor space, which is preferably the RGB color space.

The laplacian matrix L_(xy) ^(ε) follows a 2 color model specifying thatpixel color I_(i) in an input image I can be represented as a linearcombination of two colors P and S as follows: I_(i)=α_(i)P+(1−α_(i))S,∀iεw, 0≦α_(i)≦1, where the two colors are piecewise smooth and can bederived from local properties within neighborhood Ω_(x) containing thepixel i, and wherein piecewise smooth means smooth over a small region,i.e. a region similar in size as window w. It is noted that parameter εdetermines the smoothness of the decomposition into the two color modesas represented by the smoothness of α estimates. In this approach, α_(i)is defined as:

${\alpha_{i} = {{\sum\limits_{c}{a^{c}I_{i}^{c}}} + b}},{\forall{i\; \varepsilon \; w}},{c \in \left\{ {1,2,3} \right\}}$

where α are weights for each color channel, which is constant for thewindow w, b is a model parameter, and c is an index that is unque foreach window over the color channels. Preferably, the α values aredetermined according to quadratic cost function J(α)=α^(T)Lα.

In this case, the laplacian matrix L is an n×n matrix, where n is thetotal number of pixels in input image I whose (ij)^(th) element is givenby

$\sum\limits_{k|{{({i,j})} \in w_{k}}}\left( {\delta_{ij} - {\frac{1}{\left| w_{k} \right|}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{\left| w_{k} \right|}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}} \right)$

where δ_(ij) is the Kronecker delta, μ_(k) is a 3×1 mean vector ofcolors inside the k^(th) window with both i and j as members, σ_(k) isthe 3×3 covariance matrix, |w_(k)| is the cardinality of the window andI₃ is the identity matrix of size 3×3.

The laplacian matrix L may be further decomposed into a diagonal matrixD and a weight matrix W with the formulation L=D−W, wherein: diagonalmatrix D is defined with terms D_(ii)=#[k|iεw_(k)] at its diagonal andrepresents the cardinality of the number of windows of which the pixel iis a member; individual terms of weight matrix W are given by

$W_{ij} = {\sum\limits_{k|{{({i,j})} \in w_{k}}}{\frac{1}{\left| w_{k} \right|}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{\left| w_{k} \right|}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}$

where δ_(ij) is the Kronecker delta, μ_(k) is a 3×1 mean vector ofcolors inside the k^(th) window with both i and j as members, σ_(k) isthe 3×3 covariance matrix, |w_(k)| is the cardinality of the window andI₃ is the identity matrix of size 3×3; and color affinity W_(ij) for twopixels with the same color is a positive quantity varying with thehomogeneity of the local windows containing the pixels i and j, andcolor affinity for pixels with different color is zero.

Additionally at the local minima, the solution α* satisfies a firstorder optimality condition L^(T)α*=0, and the optimal condition isdefined as L^(T)α*=(D−W)^(T)α*.

The terms of W_(i,j) may be evaluated by counting the contribution ofonly the center pixel's local window defined as the window centeredabout it. Preferably, the local window w is a 3×3 pixel window requiringO(w²) computations.

In another embodiment of the present invention, the method furtherincludes steps of: (b) interpolating the input image I to a resolutionhigher than its native resolution; and (c) combining the results of step(a) with the interpolated, higher resolution image of step (b). In thisapproach, step (a) generates detail line information of the input imageI, and step (c) adds this line information to the interpolated, higherresolution image.

In the above examples, the laplacian matrix L acts like a zero-sumfilter kernel.

The present invention may further be applied to a method for enlargingan input image I. This meathod preferably includes: providing thedigital input image, I, in physical memory, the input image having colorinformation in a native color space for each image pixel; providing aprocessing device coupled to the physical memory to access the inputimage I and create an output image, h, by implementing the followingsteps: (a) applying the following sequence to the input image I,

${h_{\sigma,ɛ}(x)} = \frac{\sum\limits_{y\; \in \Omega_{x}}\left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y\; \in \Omega_{x}}\left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}$

where: (i) x and y are pixel locations over an image grid; (ii) Ω_(x) isthe neighborhood induced around the central pixel x; (iii) f_(s)(x,y)measures the spatial affinity between pixels at x and y; (iv) parameterσ controls an amount of spatial blurring; (v) L_(xy) ^(ε) is a laplacianmatrix that provides positive color affinity value for two examinedpixels, x and y, that have the same color and provides a zero coloraffinity value for pixels with different color; (vi) parameter ε is aregularization term whose relative weight determines the smoothness ofoutput image h; (vii) laplacian matrix L is defined as a diagonalmatrix, D, and a weight matrix, W, with the formulation L=D−W, wherediagonal matrix D is further defined as D_(ii)=#[k|iεw_(k)] at itsdiagonal, which represents the cardinality of the number of windows ofwhich the pixel i is a member, individual terms of the weight matrix Ware given by

$W_{ij} = {\sum\limits_{k|{{({i,j})} \in \; w_{k}}}{\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}$

and weight matrix W_(ij) is evaluated over all possible overlappingwindows that contain the center pixel; and (b) interpolating the inputimage I to a resolution higher than its native resolution; (c) combiningthe results of step (a) with the interpolated, higher resolution imageof step (b).

In this approach, step (b) generates detail line information of theinput image I, and step (c) adds this line information to theinterpolated, higher resolution image.

Preferably, step (a) is applied to the input image I in its native colorspace, which is the RGB color space.

The present invention may further be applied in a method of smoothing aninput image I, having: providing the digital input image, I, in physicalmemory, the input image having color information in a native color spacefor each image pixel; providing a processing device coupled to thephysical memory to access the input image I and create an output image,h, by implementing the following steps: (a) applying the followingsequence to the input image I,

${h_{\sigma,ɛ}(x)} = \frac{\sum\limits_{y\; \in \Omega_{x}}\left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y \in \Omega_{x}}\left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}$

where: (i) x and y are pixel locations over an image grid; (ii) Ω_(x) isthe neighborhood induced around the central pixel x; (iii) f_(s)(x,y)measures the spatial affinity between pixels at x and y; (iv) parameterσ controls an amount of spatial blurring; (v) L_(xy) ^(ε) is a laplacianmatrix that provides positive color affinity value for two examinedpixels, x and y, that have the same color and provides a zero coloraffinity value for pixels with different color; (vi) parameter ε is aregularization term whose relative weight determines the smoothness ofoutput image h; (vii) laplacian matrix L is defined as a diagonalmatrix, D, and a weight matrix, W, with the formulation L=D−W, wherediagonal matrix D is further defined as D_(ii)=#[k|iεw_(k)] at itsdiagonal, which represents the cardinality of the number of windows ofwhich the pixel i is a member, individual terms of the weight matrix Ware given by

$W_{ij} = {\sum\limits_{k|{{({i,j})} \in \; w_{k}}}{\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}$

and terms of W_(i,j) are evaluated by counting the contribution of onlythe center pixel's local window defined as the window centered about it.

Preferably, the local window w is a 3×3 pixel window. Furtherpreferably, regularization factor of ε is at last 0.1.

Like before, the present method may be applied directly to the inputimage I in its native color space, which typically is the RGB colorspace.

Other objects and attainments together with a fuller understanding ofthe invention will become apparent and appreciated by referring to thefollowing description and claims taken in conjunction with theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

In the drawings wherein like reference symbols refer to like parts.

FIG. 1 is an RGB color histogram of real world, natural, color images.

FIG. 2 illustrates a sample setup for implementing the presentBi-affinity filter.

FIG. 3 a illustrates an original image with sharp edges.

FIG. 3 b illustrates an edge-rounding effect of applying the bilateralfilter on the original image of FIG. 3 a.

FIG. 3 c illustrates an edge-preserving effect of applying theBi-affinity filter on the orignal image of FIG. 3 a.

FIG. 4 a shows a center pixel C1 and nine overlapping 3×3 pixel windowsthat include center pixel C1.

FIG. 4 b the only four overlapping windows that encompass both centerpixel C1 and target neighbor pixel N1.

FIG. 4 c shows only the center pixel's local window, which is the windowcentered about C1.

FIG. 5 illustrates the affect of varying regularization term ε pixelwindows of varying pixel size.

FIG. 6 shows the effect of regularization term ε for a fixed windowsize, as regularization term ε assigned values of 0.0005, 0.005, 0.05,and 0.5.

FIGS. 7 a and 7 b show comparisons of the PSNR of the presentBi-affinity filter applied in an image's native RGB color space versusthe bilateral filter applied in both the image's native RGB color spaceand converted CIELAB color space are shown for a fixed range filtervariance and a varied window size.

FIG. 8 shows compares the results of applying the present Bi-affinityfilter to an input image with the results of applying a bilatera filterto the input image in its native RGB color domain and after beingconverted to the CIELAB color space.

FIG. 9 shows an additional set of experimental results similar to thoseof FIG. 8 for comparing the Bi-affinity filter to the bilateral filterin both the native RGB color space and in the CIELAB color space.

FIG. 10 shows an additional set of experimental results similar to thoseof FIG. 8 for comparing the Bi-affinity filter to the bilateral filterin both the native RGB color space and in the CIELAB color space.

FIGS. 11 a, 11 b, and 11 c show the application of the presentBi-affinity filter to enhance a low resolution image.

FIG. 12 a shows an image zoomed-in by a factor of 2× by bi-cubicinterpolation.

FIG. 12 b shows the samed zoomed-in image of FIG. 12 a additionallyenhanced by the present Bi-affinity filter.

FIG. 13 a shows an image zoomed-in by a factor of 2× by bi-cubicinterpolation.

FIG. 13 b shows the zoomed-in image of FIG. 13 a additionally enhancedby the present Bi-affinity filter.

FIG. 14 a shows an image zoomed-in by a factor of 2× by bi-cubicinterpolation.

FIG. 14 b shows the zoomed-in image of FIG. 14 a additionally enhancedby the present Bi-affinity filter.

FIGS. 15 a and 15 b show a graphical model for image enhancement.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

It is presently proposed that the true power of the bilateral filter isyet to be realized owing to the fact that very few family of kernelshave been tried with the bilateral paradigm.

Traditionally, research has overlooked one of the most importantshortfalls of the bilateral filter, which is a unified handling ofmulti-channel color images. This is due to the fact that color channelsare traditionally assumed to be independent of each other, and thefilter therefore processes each color channel independently. Forexample, an RGB color space, which has three color channels (red, green,and blue), requires separate application of the bilateral filter to eachof three color channels. As a direct consequence of this, the bilateralfilter produces color artifacts at sharp color edges in the RGB colorspace. This poses a problem since the RGB color space is typically thenative color space of color digital images, i.e. images obtained throughvarious digital imaging techniques, such as digital photography or imagescanning operations.

One remedy previously proposed is to convert from the native RGB colorspace to the CIELAB color space, which attempts to approximate humanvision (or perception). Like the RGB color space, the CIELAB color spacealso has three color channels, and thus still requires that thebilateral filter be applied independently to each of three colorchannels, but the CIELAB is noted for avoiding color artifacts. That is,according to this approach, once the native RGB color image is convertedto the CIELAB color space, then channel wise bilateral filter processingdoes not produce the artifacts apparent from processing in the RGB colorspace.

The present discussion further investigates this weakness and proposes anew technique that works at par with the transformed domain techniques(i.e. converting from a native RGB color space to an alternate colorspace) that have been the standard practices within the digital imagingcommunity, thus far.

A new filter, hereafter called a “Bi-affinity filter”, for color imagesis presently proposed. This filter is similar in structure to abilateral filter, but the proposed Bi-affinity filter is based on acolor line model. As will be shown, this approach permits theelimination of the explicit conversion from the RGB color space toperception based color spaces, such as CIELAB.

The present Bi-affinity filter measures the color affinity of a pixel toa small neighborhood around it and weighs the filter term accordingly.It is put forth that this method can perform at par with standardbilateral filters for color images. An added benefit of the presentBi-affinity filter is that not only does it preserve big edges of aninput image (a desired characteristic of the bilateral filter), but alsopreserves and enhances small edges of the input image, which leads toadditional applications as an image enhancement filter.

As is explained above, the principle objective of the bilateral filteris to combine information from the spatial domain with a feature domain,and this feature domain is generally the color domain. Determining coloraffinity between pixels in color images in the RGB domain is complicatedby differences in light intensity causing pixels that are close in colorto appear quite different and thus have markedly different RGB values.It is therefore generally preferred to convert from the RGB color spaceto an alterante color space better capable to filtering out differencein light intensity, and apply the bilateral filter in the alternatecolor space.

Many alternative color spaces have been suggested to separate color fromlight intensity and create methods for determining whether two pixelsshare the same color, or not. These color spaces can be divided into twogroups, namely linear and non-linear. Among the linear color spaces, themost widely used are the YCrCb, YUV, and YIQ color spaces. Among thenon-linear color spaces, two sub-groups are popular: a first sub-groupthat separates color into hue (i.e. color), saturation (i.e. purity),and value (i.e. light intensity), and a second sub-group that separatescolor into luminance plus two color coordinates. The first group mayfurther include the HSV, HSI, HSB, HSL, etc. color spaces. The secondgroup are the CIELAB and CIE-LUV color spaces that separate color intoluminance (i.e. light intensity) and two color coordinates in an effortto create a color space that is perceptually uniform (i.e. close to howthe human eye perceives color variations).

To determine whether two pixels have the same real world color, or not,the color coordinates of a generic color model are used. These colormodels assume either that there is no color distortion or that there isan identical color distortion for all imaging conditions. In practice,when dealing with real world images of an unknown source, theseassumptions are rarely true as scene surface color is distorteddifferently in different images as well as different image regions,depending on the scene and camera settings.

This leads to the topic of the color line model. The introduction ofcolor lines has been attributed to Omer et al. in Color lines: ImageSpecific Color Representation, CVPR, 2004, herein incorporated in itsentirety by reference. Omer et al. proposed that color natural imagesform clusters of pixel colors in the RGB color space. These clusters ofpixels appear to form mostly tubular regions, due to the fact that mostsmall regions in natural images can be decomposed into a linearcombination of 2 colors. This lead to the two-color model, describedbelow.

With reference to FIG. 1, when looking at the RGB color histogram 10 ofreal world, natural images, it can be clearly observed that thehistogram is very sparse and structured. The color line model exploitsthese two properties of RGB color histograms by describing the elongatedcolor clusters, or tubular regions, 11-15. Since the structure andarrangment of the color clusters 11-15 are specific to an image, thisresults in an image-specific color representation that has two importantproperties: robustness to color distortion and a compact description ofcolors in an image. This idea has been used for image matting (i.e.separation of a foreground subject from its surrounding background),Bayer demosaicing (i.e. reconstruction of a color image from incompletecolor samples output from an image sensor overlaid with a Bayer colorfilter array), and more recently for image de-noising and de-blurring.

The two-color model states that any pixel color I_(i) in an image I canbe represented as a linear combination of two colors P and S, wherethese colors are piecewise smooth and can be derived from localproperties within a small neighborhood containing the pixel i.

I _(i)=α_(i) P+(1−α_(i))S, ∀iεw, 0≦α_(i)≦1  (5)

where w is a small patch, or window, and α is the matting coefficient.The patch size is a key parameter in this two-color model, as it is truefor only small neighborhoods, relative to the resolution and size of theimage. These “small” neighborhoods are typically defined as 3×3 pixelneighborhoods. As the resolution and the size of the images grow, soshould the window size as well, to capture a valid neighborhood.

As Levin et al. show in A Closed Form Solution to Natural Image Matting,CVPR, 2006 (herein incorporated in its entirty by reference), a similarline of reasoning may be used for color images. For color images, if thecolor line property is obeyed, then the 4D linear model satisfied by thematting coefficient, within a small window, at each pixel can be writtenas

$\begin{matrix}{{\alpha_{i} = {{\sum\limits_{c}{a^{c}I_{i}^{c}}} + b}},{\forall\; {i \in w}},{c \in \left\{ {1,2,3} \right\}}} & (6)\end{matrix}$

Where α are weights for each color channel, which is constant for thewindow w, b is a model parameter, and c is an index that is unque foreach window over the color channels.

A cost function for evaluating the matting coefficient α can beformulated from this model. For an image with N pixels, the cost isdefined as

$\begin{matrix}{{J\left( {\alpha,a,b} \right)} = {\sum\limits_{k \in N}\left( {{\sum\limits_{i \in w_{k}}\left( {\alpha_{i} - {\sum\limits_{c}{a_{k}^{c}I_{i}^{c}}} - b_{k}} \right)^{2}} + {ɛ{\sum\limits_{c}a_{k}^{c^{2}}}}} \right)}} & (7)\end{matrix}$

where w_(k) is a small window around pixel k and a={α_(i) ^(c)}, for alli=[1,N]. ε is a regularization weight for uniqueness as well assmoothness of the solution.

A first theorem is herein proposed, as follow.

Theorem 1: Let J(α)=min_(a,b)J(α,a, b), then J(α)=α^(T)Lα, where L is alaplacian matrix. As it is known in the art, and in particular in themathematical field of graph theory, a laplacian matrix (also called anadmittance matrix or Kirchhoff matrix, is a matrix representation of agraph. Together with Kirchhoff's theorem, it can be used to calculatethe number of spanning trees for a given graph. In the present case, Lis an n×n matrix, whose ij^(th) element is given by

$\begin{matrix}{\sum\limits_{k|{{({i,j})} \in \; w_{k}}}\left( {\delta_{ij} - {\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)^{T}{{\overset{\sim}{\Sigma}}_{k}^{- 1}\left( {I_{j} - \mu_{k}} \right)}}} \right)}} \right)} & (8)\end{matrix}$

where δ_(ij) is the Kronecker delta, μ_(k) is a 3×1 mean vector ofcolors inside the k^(th) window with both i and j as members, I_(i) andI_(j) are the color vectors at location i and j,

${{\overset{\sim}{\Sigma}}_{k} = {\Sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}}},$

where Σ_(k) is the 3×3 covariance matrix, |w_(k)| is the cardinality ofthe window and I₃ is the 3×3 identity matrix.Proof: Levin et. al. provide a limited proof based on an extension froma gray scale case. Below is presented a full 3-channel proof that can bereadily extended to more channels, if necessary. Rewriting Eq. 7 in amatrix notation, where ∥•∥ denotes the 2-norm,

$\begin{matrix}{{{J\left( {\alpha,a,b} \right)} = {\sum\limits_{k}{{{G_{k} \cdot {\overset{\_}{a}}_{k}} - \alpha_{k}}}}}{where}} & (9) \\{{G_{k} = \begin{bmatrix}I_{1}^{R} & I_{1}^{G} & I_{1}^{B} & 1 \\\vdots & \vdots & \vdots & \vdots \\I_{w_{k}}^{R} & I_{w_{k}}^{G} & I_{w_{k}}^{B} & 1 \\\sqrt{ɛ} & 0 & 0 & 0 \\0 & \sqrt{ɛ} & 0 & 0 \\0 & 0 & \sqrt{ɛ} & 0\end{bmatrix}},{{\overset{\_}{a}}_{k} = \begin{bmatrix}a_{k}^{R} \\a_{k}^{G} \\a_{k}^{B} \\b\end{bmatrix}},{\alpha_{k} = \begin{bmatrix}\alpha_{1} \\\vdots \\\alpha_{w_{k}} \\0 \\0 \\0\end{bmatrix}}} & (10)\end{matrix}$

Note that another representation of G_(k) is possible where the last 3rows are combined to a single row of the form [√{square root over(ε)}√{square root over (ε)}√{square root over (ε)}0], but this formleads to an unstable covariance matrix. For known α_(k), the leastsquare problem can be solved

ā _(k)=argmin∥G _(k) ·ā _(k)−α_(k)∥  (11)

=(G _(k) ^(T) G _(k))⁻¹ G _(k) ^(T)α_(k)  (12)

Substituting this solution in Eq. 9, and denoting L_(k)=I_(|w) _(k)_(|+3)−G_(k)(G_(k) ^(T)G_(k))⁻¹G_(k) ^(T), where I_(|w) _(k) _(|+3) isthe identity matrix of size (|w_(k)|+3), one obtains,

$\begin{matrix}\begin{matrix}{{J(\alpha)} = {\sum\limits_{k}{{L_{k}\alpha_{k}}}}} \\{= {\sum\limits_{k}\left( {\alpha_{k}^{T}L_{k}^{T}L_{k}\alpha_{k}} \right)}}\end{matrix} & (13)\end{matrix}$

Making the additional observation that

$\begin{matrix}{{L_{k}^{T}L_{k}} = {\left( {I_{{w_{k}} + 3} - {{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}G_{k}^{T}}} \right)^{T}\left( {I_{{w_{k}} + 3} - {{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}G_{k}^{T}}} \right)}} \\{= {I_{{w_{k}} + 3} + {{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}G_{k}^{T}{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}G_{k}^{T}} - {2{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}G_{k}^{T}}}} \\{= {I_{{w_{k}} + 3} - {{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}G_{k}^{T}}}} \\{= L_{k}}\end{matrix}$

one can write equation (13) as J(α)=Σ_(k)(α_(k) ^(T)L_(k)α_(k)). Tocomplete the proof one needs to find the expression for L_(k)|_(i,j).

Noting the identity E[X²]=σ_(XX) ²+E[X]², denoting the individualchannel means E[R] as R, one can write

$\begin{matrix}{{G_{k}^{T}G_{k}} = \mspace{675mu} (14)} \\{\; {{w}\begin{bmatrix}\overset{\overset{A}{}}{\left( \begin{matrix}\begin{matrix}\begin{matrix}{\sigma_{RR}^{2} + R^{2} + \frac{ɛ}{w_{k}}} \\{\sigma_{GR}^{2} + {GR}}\end{matrix} \\{\sigma_{BR}^{2} + {BR}}\end{matrix} & \begin{matrix}\begin{matrix}{\sigma_{RG}^{2} + {RG}} \\{\sigma_{GG}^{2} + G^{2} + \frac{ɛ}{w_{k}}}\end{matrix} \\{\sigma_{BG}^{2} + {BG}}\end{matrix} & \left. \begin{matrix}\begin{matrix}{\sigma_{RB}^{2} + {RB}} \\{\sigma_{GB}^{2} + {GB}}\end{matrix} \\{\sigma_{BB}^{2} + B^{2} + \frac{ɛ}{w_{k}}}\end{matrix} \right)\end{matrix} \right.} & \overset{\overset{D}{}}{\begin{pmatrix}\begin{matrix}\underset{\;}{R} \\\underset{\;}{\overset{\;}{G}}\end{matrix} \\\overset{\;}{B}\end{pmatrix}} \\\underset{\underset{D^{T}}{}}{\begin{matrix}\begin{matrix}\left( R \right. & \; & \; & \; & \mspace{11mu} & \; & \; & \; & \; & \; & G\end{matrix} & \; & \; & \; & \; & \mspace{11mu} & \mspace{11mu} & \; & \; & \; & \left. B \right)\end{matrix}} & \underset{\underset{C}{}}{\overset{\;}{1}}\end{bmatrix}}}\end{matrix}$

where the matrix has been divided into 4 components. Note that D=μ_(k)for the k^(th) window. The inverse of the above system can now bewritten as:

$\left( {G_{k}^{T}G_{k}} \right)^{- 1} = {\frac{1}{w_{k}}\begin{bmatrix}P & Q \\R & S\end{bmatrix}}$ $\begin{matrix}{P = \left( {A - {{DC}^{- 1}D^{T}}} \right)^{- 1}} \\{= \left( {A - {DD}^{T}} \right)^{- 1}} \\{= \begin{bmatrix}{\sigma_{RR}^{2} + \frac{ɛ}{w_{k}}} & \sigma_{RG}^{2} & \sigma_{RB}^{2} \\\sigma_{GR}^{2} & {\sigma_{GG}^{2} + \frac{ɛ}{w_{k}}} & \sigma_{GB}^{2} \\\sigma_{BR}^{2} & \sigma_{BG}^{2} & {\sigma_{BB}^{2} + \frac{ɛ}{w_{k}}}\end{bmatrix}^{- 1}} \\{= {\overset{\sim}{\Sigma}}_{k}^{- 1}}\end{matrix}$$Q = {{- {P\left( {DC}^{- 1} \right)}} = {{- {PD}} = {{- {\overset{\sim}{\Sigma}}_{k}^{- 1}}\mu_{k}}}}$$R = {{{- \left( {C^{- 1}D^{T}} \right)}P} = {{{- D^{T}}P} = {{- \mu_{k}^{T}}{\overset{\sim}{\Sigma}}_{k}^{- 1}}}}$$S = {{C^{- 1} - {R\left( {DC}^{- 1} \right)}} = {{1 - {RD}} = {1 + {\mu_{k}^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}\mu_{k}}}}}$

Putting all the terms together, one can write

$\begin{matrix}{\left( {G_{k}^{T}G_{k}} \right)^{- 1} = {\frac{1}{w_{k}}\begin{bmatrix}{\overset{\sim}{\Sigma}}_{k}^{- 1} & {{- {\overset{\sim}{\Sigma}}_{k}^{- 1}}\mu_{k}} \\{{- \mu_{k}^{T}}{\overset{\sim}{\Sigma}}_{k}^{- 1}} & {1 + {\mu_{k}^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}\mu_{k}}}\end{bmatrix}}} & (15) \\{{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1} = {\frac{1}{w_{k}}\begin{bmatrix}{\left( {I_{1} - \mu_{k}} \right)^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}} & {1 - {\left( {I_{1} - \mu_{k}} \right)^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}\mu_{k}}} \\{\left( {I_{2} - \mu_{k}} \right)^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}} & {1 - {\left( {I_{2} - \mu_{k}} \right)^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}\mu_{k}}} \\\vdots & \vdots \\{\left( {I_{w_{k}} - \mu_{k}} \right)^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}} & {1 - {\left( {I_{w_{k}} - \mu_{k}} \right)^{T}{\overset{\sim}{\Sigma}}_{k}^{- 1}\mu_{k}}} \\{\sqrt{ɛ}{\overset{\sim}{\Sigma}}_{k}^{- 1}} & {\sqrt{ɛ}{\overset{\sim}{\Sigma}}_{k}^{- 1}\mu_{k}}\end{bmatrix}}} & (16)\end{matrix}$

Right multiplication by G_(k) ^(T) yields the final symmetric form,where only the i^(th) column is shown for conciseness and ease ofunderstanding

${{G_{k}\left( {G_{k}^{T}G_{k}} \right)}^{- 1}{G_{k}^{T}\left\lbrack {:{,i}} \right\rbrack}} = {\frac{1}{w_{k}}\begin{bmatrix}{1 + {\left( {I_{1} - \mu_{k}} \right)^{T}{{\overset{\sim}{\Sigma}}_{k}^{- 1}\left( {I_{i} - \mu_{k}} \right)}}} \\{1 + {\left( {I_{2} - \mu_{k}} \right)^{T}{{\overset{\sim}{\Sigma}}_{k}^{- 1}\left( {I_{i} - \mu_{k}} \right)}}} \\{1 + {\left( {I_{3} - \mu_{k}} \right)^{T}{{\overset{\sim}{\Sigma}}_{k}^{- 1}\left( {I_{i} - \mu_{k}} \right)}}} \\\vdots \\{1 + {\left( {I_{w_{k}} - \mu_{k}} \right)^{T}{{\overset{\sim}{\Sigma}}_{k}^{- 1}\left( {I_{i} - \mu_{k}} \right)}}} \\{ɛ{{\overset{\sim}{\Sigma}}_{k}^{- 1}\left( {I_{i} - \mu_{k}} \right)}}\end{bmatrix}}$

Subtracting from I_(|w) _(k) _(|+3) and summing over k concludes theproof. Note that G_(k) has 3 extra rows (or C extra rows for generalcase) for the regularization ε. These can be neglected in the finalexpression since they do not explicitly effect the other computations.

As is stated above, laplacian matrix L, is an n×n matrix. It is notethat in L, n may also be understood to be the total number of pixels inthe image, whose (ij)^(th) element is given by

$\begin{matrix}{\sum\limits_{k{{({i,j})} \in w_{k}}}\; \left( {\delta_{ij} - {\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}} \right)} & (17)\end{matrix}$

where δ_(ij) is the Kronecker delta, μ_(k) is a 3×1 mean vector ofcolors inside the k^(th) window with both i and j as members, σ_(k) isthe 3×3 covariance matrix, |w_(k)| is the cardinality of the window andI₃ is the identity matrix of size 3×3. Note that the term ε is aregularization term and determines the smoothness of the decompositioninto two dominant color modes (i.e. two colors P and S). The laplacianmatrix defined by L has been termed the matting laplacian.

The usual decomposition of the laplacian matrix L into a diagonalmatrix, D, and a weight matrix, W, leads to the formulation L=D−W. HereD is a diagonal matrix with the terms D_(ii)=#[k|iεw_(k)] at itsdiagonal, which represents the cardinality of the number of windows ofwhich the pixel i is a member. The individual terms of the weight matrixW are given by

$\begin{matrix}{W_{ij} = {\sum\limits_{k{{({i,j})} \in w_{k}}}{\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}} & (18)\end{matrix}$

This brings the present discussion to the subject of the Bi-affinityfilter. It is first noted that by the definition of the laplacianmatrix, all its rows sum to zero, which leads to D_(ii)=Σ_(j)W_(ij). Atthe local minima, the solution α* satisfies the first order optimalitycondition L^(T)α*=0. So the optimal condition can be written as

$\begin{matrix}{{L^{T}\alpha^{*}} = {\left( {D - W} \right)^{T}\alpha^{*}}} & (19) \\{{L^{T}\alpha^{*}} = \begin{pmatrix}{D_{11}\alpha_{1}^{*}} & {- {\sum\limits_{j}\; {W_{1\; j}\alpha_{j}^{*}}}} \\{D_{22}\alpha_{2}^{*}} & {- {\sum\limits_{j}\; {W_{2\; j}\alpha_{j}^{*}}}} \\\vdots & \vdots \\{D_{nn}\alpha_{n}^{*}} & {- {\sum\limits_{j}\; {W_{Nj}\alpha_{j}^{*}}}}\end{pmatrix}} & (20)\end{matrix}$

Substituting D_(ii)=Σ_(j)W_(ij) into the above system of equations andinvoking the first order optimality condition leads to

$\begin{matrix}{\begin{pmatrix}{\sum\limits_{j}\; {\left( {\alpha_{1}^{*} - \alpha_{j}^{*}} \right)W_{1\; j}}} \\{\sum\limits_{j}\; {\left( {\alpha_{2}^{*} - \alpha_{j}^{*}} \right)W_{2\; j}}} \\\vdots \\{\sum\limits_{j}\; {\left( {\alpha_{n}^{*} - \alpha_{j}^{*}} \right)W_{nj}}}\end{pmatrix} = 0} & (21)\end{matrix}$

The effect of this equation is that the color affinity W_(ij) for twopixels with the same color (same α*), is a positive quantity varyingwith the homogeneity of the local windows containing the pixels i and jas governed by Eqn. 18. But for pixels with different color (differentα*) the color affinity is zero. In essence the rows of the laplacianmatrix L work as a zero-sum filter kernel, after appropriate resizing ofthe window.

This leads to the formulation of the Bi-affinity filter as:

$\begin{matrix}{{h_{\sigma,ɛ}(x)} = \frac{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}} & (22)\end{matrix}$

FIG. 2 illustrates a sample setup for implementing the presentBi-affinity filter. The input image may be provided in a digital memory110 that is in communication with a data processing device 112, such asa CPU, ASIC, PLD, CPLD, or other type of data processing or computingdevice. In the present case, the data processing device processing theinput image in memory 110 according to processing steps defined byequation 22, where the dependence on user specified parameters σ,ε onthe filter output are denoted. I and h are the input and output imagesrespectively, x and y are pixel locations over an image grid, and Ω_(x)is the neighborhood induced around the central pixel x. In the presentcase, f_(s)(x,y) measures the spatial affinity between pixels at x andy, where parameter σ controls the amount of spatial blurring and servesa similar purpose as the spatial filter variance (i.e. σ_(s) in eqn. 3)in a traditional bilateral filter. L_(xy) ^(ε) is the laplacian matrix(or matting laplacian) that provides positive color affinity value fortwo examined pixels, x and y, that have the same color and provides azero color affinity value for pixels with different color. The parameterε, a regularization term, works analogous to the range varianceparameter (i.e. σ_(r) in eqn. 4) in a traditional bilateral filter. Notethat the relative weight attributed to regularization term ε, determinesthe smoothness of the α estimates (i.e. the smoothness of color blendingfrom pixel to pixel), which in the present work translates to thesmoothness of the filtered image. The output image h may be saved backto memory 110 or may be stored in another memory device, not shown.

However, the Bi-affinity filter does not smooth the edge, due to theaffinity formulation which is zero across the edge. This is in directcontrast to the Bilateral filter, which has an inherent bending effectat the edges. This image distortion can be observed in FIGS. 3 a-3 c.FIG. 3 a is an example of an origial image with a sharp image. FIG. 3 bshows the result of applying bilateral filtering. As is shown, bilateralfiltering results in rounding of the edges. By contrast, applyingBi-affinity filtering preserves edge sharpness, as is shown in FIG. 3 c.

The calculation of the exact weight matrix W_(ij) (hereafter called the“affinity matrix”) as mentioned in Eqn. 18, involves evaluations overall possible overlapping windows that contain the center pixel. Thisinvolves O(w³) computations, where w is the size of the window. Forexample, if the window size is 3 (i.e. three pixels by three pixels),then nine computations would be required per center pixel, which is morethan typically required in a traditional bilatera filter.

For example FIG. 4 a shows a center pixel C1 and nine overlapping 3×3pixel windows 21-29 that include center pixel C1. Neighbor pixels N1-N8are also noted. In the present example, the color affinity betweencenter pixel C1 and target neigbhor pixel N1 is desired.

The overall complexity can be reduced by evaluating the color affinityover a smaller set of possible windows. For example in FIG. 4 b, onlythe four overlapping windows, 22, 23, 24, and 29, that encompass bothcenter pixel C1 and target neighbor pixel N1 may be evaluated.

In the simplest case, shown in FIG. 4 c, the terms of W_(i,j) can beevaluated locally, i.e. counting the contribution of only the centerpixel's local window, 29, which is the window centered about C1. In thiscase, the complexity is similar to the typical bilateral range filter,which involves O(w²) compuations. This simplest case of only evaluatingterms of W_(i,j) local to a center pixel (i.e. only evaluating thecenter pixel's local window) is the basis of an approximate Bi-affinityfilter.

To keep fair the following comparisons of the Bi-affinity filter withthe bilateral filter (i.e. to assure that both both have a similar levelof complexity), the bilatera filter will instead be compared to theapproximate Bi-lateral filter denoted by h^(l)(x) and herein defined as:

$\begin{matrix}{{h^{l}(x)} = \frac{\sum\limits_{y \in \Omega_{x}}\; {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{x\; ɛ}{I(y)}}}{\sum\limits_{y \in \Omega_{x}}\; {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{x\; ɛ}}}} & (23)\end{matrix}$

which considers only the local window centered around pixel x, denotedby L^(x).

The operations involved in computing the L_(ij) terms as mentioned inEqn. 6 (or Eqn. 17), can be decomposed as summation of Gaussianlikelihoods over window dependent parameters μ_(w), Σ_(w). Theseparameters can be computed by accumulating first and second ordersufficient statistics over windows. If memory complexity is not an issuethen pre-computing 9 integral images can be an option. These 9 integralimages correspond to 3 integral images for each of the R, G and B colorchannels, 3 integral images for RR, GG and BB channels, and 3 integralimages for RG, GB and RB channels. For 3 channel color images, this isequivalent to storing 3 more images into the memory. For really largeimages (HDTV etc.) this option might not be the most optimal due to thelarge memory overhead.

Another method of computing the L_(ij) terms is to collect sufficientstatistics for the current window and then update the statistics foreach unit move from top to bottom and left to right, as proposed by themedian filtering approach in Transforms and Median Filters (inTwo-Dimensional Signal Processing II, pages 209-211, Berlin, 1981.Springer-Verlag, Huan), which is herein incorporated in its entirety byreference and in Fast Median and Bilateral Filtering (ACM Trans. Graph.,25(3):519-526, 2006, by Weiss), which is also herein incorporated in itsentirety by reference. Both of these methods can now be used toimplement both the full and the approximate bi-affinity filter.

The above approximate Bi-affinity filter (as described in eqn. 23) wastested.

With reference to FIG. 5, the regularization term ε in the affinityformulation works as an edge smoothness term. For understanding theeffect of this term, the amount of regularization ε used for the processwas varied for each of nine square windows, W5, W7, W9, W11, W13, W15,W17, W19, and W21, with each window having a size of 5, 7, 9, 11, 13,15, 17, 19, and 21 pixels per side, respectively. The power-to-signalnoise ratio, PSNR, for each window with respect to the original imagewas recorded. As shown, the PSNR degrades for larger window sizes. Thisis in keeping with the two color model, which is valid only for smallwindows. However, the regularization term ε neutralizes the effect ofwindow size to a certain degree, as is evident by the band of valuescollecting near the PSNR value of 96 dB. This suggests a possibletradeoff between PSNR and edge smoothness. For very small regularizationvalues, the noise across the edge can contribute to the jaggedness ofthe reconstruction. This effect can be countered by increasingregularization term ε, but the increased regularization comes at thecost of increased smoothness of the overall image.

Empirically, good results have been obtained for larger window sizes bykeeping values of regularization terms relatively larger than thoseproposed in the matting literature mentioned above (i.e. preferaby atleast 1.5 times larger).

With reference to FIG. 6, the effect of regularization term ε for afixed window size are illustrated by showing the results assigningvalues 0.0005, 0.005, 0.05, and 0.5 to regularization term ε. The edgereconstruction becomes increasingly jagged as the amount ofregularization is decreased, as is more clealry seen in the closeup viewat the upper-right corner of each sample image.

With reference to FIGS. 7 a and 14 b, comparisons of the PSNR of thepresent Bi-affinity filter applied in an image's native RGB color spaceversus the bilateral filter applied in both an image's native RGB colorspace and converted CIELAB color space are shown. For quantitativecomparisons against the traditional bilateral filter, the window size isvaried for a fixed range filter variance of 0.1 in FIG. 7 a, and for afixed range filter variance of 1.0 in FIG. 7 b.

With ε=σ_(r)=0.1 in FIG. 7 a, the PSNR values of the present Bi-affinityfilter 100 is better than that of the RGB bilateral filter 102, but islower than the CIELAB bilateral filter 104. Nonetheless, the presentBi-affinity filter 100 is still within acceptable deviatipons from theCIELAB bilatera filter 104.

However, as is shown in FIG. 7 b, when the range filter variance isincreased to 1.0 (i.e. ε=1), the perfance of the present Bi-affinityfilter 100 can surpass that of the CIELAB bilatera filter 104. Morespecifically, when the window is increase to at least 5×5 pixels (i.e.w=5), the performance of the present Bi-affinity filter in the image'snative RGB color spaces surpasses that of the bilateral filter even inthe converted CIELAB color space.

As is mentioned earlier in reference to Eqn. 23, the approximateBi-affinity filter, which is evaluated over one centered local window,closely emulates the traditional bilateral filter. The following, nextexperimental results compare the present approximate Bi-affinity filterto the traditional bilateral filter. These results are illustrated inFIG. 8.

With reference to FIG. 8, the upper-left corner image shows the originalinput image, to which the bilateral filter and the present approximteBi-affinity filter will be applied. The results of applying thebileteral filter in the input image's native RGB color space (with thebilateral filter individually applied to each of the Red color channel,Green color channel, and Blue color channel) is shown in the upper-rightcorner. As shown, the resultant soften image suffers from colorartifacts at sharp color edges that may be perceived as a loss in imagedefinition.

However, if the input image is first converted from the RGB color spaceto the CIELAB color space, and the bilateral filter is then individuallyapplied to each of the three CIELAB color channels, then the resultantsoften image is more clearly defined, as is shown in the lower-leftcorner of FIG. 8.

Nonetheless, as is shown in the lower-right corner of FIG. 8, similarresults are obtained by applying the present Bi-affinity filter to theoriginal input image in its native RGB color space, with aregularization factor of ε=0.1. Thus, the response of the bilateralfilter in the CIELAB color space and the Bi-affinity filter in thenative RGB color space is very similar, even though the Bi-affinityfilter does not need, or use, any color conversions, resulting in areduction in filter processing steps.

Two additional sets of experimental results similar to those of FIG. 8for comparing the Bi-affinity filter to the bilateral filter in both thenative RGB color space and in the CIELAB color space are shown in FIGS.9 and 10.

The present Bi-affinity filter also has applications in enhanced imagezooming. The original Bi-affinity filter of Eqn. 22 (i.e. the versionnot limited to evaluation over only one centered local window), has beenshown to preserve very minute details. In other words, the presentbi-affinity formulation preserves very intricate details of an image ascompared to the traditional bilateral filter, which preserves onlydominant edges of an image. In this regard, the present Bi-affinityfilter can be understood to preserve all edges, whereas the bilateralfilter preserves only strong edges. This important feature of theBi-affinity filter leads to one of the most interesting applications ofthe Bi-affinity filter: image enhancement and zooming.

Image enhancement techniques attempt to convert a low-resolution imageto a high-resolution image. In effect, image enhancment techniques tryto generate high-resolution data from available low-resolution data byestimating missing information. This leads to numerous formulations,some learning based and some interpolation based. Basically, if one caninfer the missing high-resolution data, then one can add the missinghigh-resolution data to an interpolation of the low-resolution data(i.e. an enlarged or zoomed-in version of the input image) thatsatisfies data fidelity constraints. The result is an estimated (i.e.generated) high-resolution image.

For example with reference to FIG. 11 a, given a low-resolution inputimage, as shown, one can interpolate it to a desired higher-resolution,and then add the missing high-resolution information to generate a finalhigh-resolution result. The mean affinity at each pixel, which is therow wise normalized summation of W, appears to contain this missingdetail. This missing detail is shown in FIG. 11 b. The Bi-affinityfilter further places a smoothed affinity weighted kernel at each pixel,and the combined effect is the enhanced image shown in FIG. 11 c. Thisenhancement effect is a bi-product of the filtering formulation of theBi-affinity filter.

A second example of the enhancement feature of the present Bi-affinityfilter is shown FIGS. 12 a and 12 b. FIG. 12 a shows an image zoomed-inby a factor of 2× by bi-cubic interpolation. FIG. 12 b shows the samedzoomed-in image additionally enhanced by the present Bi-affinity filter.Note the preservation of small image details.

A third example of the enhancement feature of the present Bi-affinityfilter is shown FIGS. 13 a and 13 b. FIG. 13 a shows an image zoomed-inby a factor of 2× by bi-cubic interpolation. FIG. 13 b shows the samedzoomed-in image additionally enhanced by the present Bi-affinity filter.

A fourth example of the enhancement feature of the present Bi-affinityfilter is shown FIGS. 14 a and 14 b. FIG. 14 a shows an image zoomed-inby a factor of 2× by bi-cubic interpolation. FIG. 14 b shows the samedzoomed-in image additionally enhanced by the present Bi-affinity filter.

A graphical model for image enhancement is shown in FIGS. 15 a and 15 b.FIG. 15 a illustrates a passive filtering technique and FIG. 15 billustrates an active filtering technique. In both FIGS. 15 a and 15 b,the circles represent observed nodes, the squares represent functionnodes, and the triangles represent hidden nodes to be estimated.

The method of the present invention is a passive filtering technique,and is best described by FIG. 15 a. Like other passive filteringtechniques, e.g. bilateral, bicubic, etc., the present method only looksat the low-resolution observation layer (i.e. layer having lower circlesLC along the dashed grid pattern) to generate the values of thehigh-resolution scene (i.e. the upper layer having upper triangles UT).In the present case, the filter kernel is applied to the observed nodeto generate the scene.

By comparison, active methods such as Markov random field (MRF) basedmodels, shown in FIG. 15 b, impose neighborhood continuity constraintsin the inferred layer as well. That is, scene nodes are also influencedby the neighbors in the hidden scene field. This is typically donethrough interaction potential within the scene nodes. It is presentlyput forth that the principles of the present invention may be applied toa such model, as well.

In summary, a new edge preserving filter, which works on the principleof matting affinity, is proposed. The formulation of matting affinityallows a better representation of the range filter term in bilateralfilter class. The definition of the affinity term can be relaxed to suitdifferent applications.

An approximate Bi-affinity filter whose output is shown to be verysimilar to the traditional bilateral filter is defined. The presenttechnique has the added advantage that no color space changes arerequired and hence an input images can be handled in its original colorspace, i.e. its native color space. This is a big benefit overtraditional bilateral filter, which needs a conversion from its nativecolor space to a perception based color spaces, such as CIELAB, togenerate results close to the present invention.

Furthermore, the full Bi-affinity filter preserves very minute detailsof the input image, and can easily be extended to an image enhancementapplication to permit an enhanced zooming function.

While the invention has been described in conjunction with severalspecific embodiments, it is evident to those skilled in the art thatmany further alternatives, modifications, and variations will beapparent in light of the foregoing description. Thus, the inventiondescribed herein is intended to embrace all such alternatives,modifications, applications and variations as may fall within the spiritand scope of the appended claims.

1. A method of processing a digital, input image, comprising: providingsaid digital input image, I, in physical memory, said input image havingcolor information in a native color space for each image pixel;providing a processing device coupled to said physical memory to accesssaid input image I and create an output image, h, by implementing thefollowing steps: (a) applying the following sequence to said input imageI,${h_{\sigma,ɛ}(x)} = \frac{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}$where x and y are pixel locations over an image grid, and Ω_(x) is theneighborhood induced around the central pixel x; f_(s)(x,y) measures thespatial affinity between pixels at x and y; parameter σ controls anamount of spatial blurring; L_(xy) ^(ε) is a laplacian matrix thatprovides positive color affinity value for two examined pixels, x and y,that have the same color and provides a zero color affinity value forpixels with different color; parameter ε is a regularization term whoserelative weight determines the smoothness of output image h.
 2. Themethod of claim 1, wherein step (a) is applied to said input image I inits native color space.
 3. The method of claim 2, wherein said nativecolor space is the RGB color space.
 4. The method of claim 1, whereinsaid laplacian matrix L_(xy) ^(ε) follows a 2 color model specifyingthat pixel color I_(i) in an input image I can be represented as alinear combination of two colors P and S as follows:I _(i)=α_(i) P+(1−α_(i))S, ∀iεw, 0≦α_(i)≦1 where the two colors arepiecewise smooth and can be derived from local properties withinneighborhood Ω_(x) containing the pixel i, and wherein piecewise smoothmeans smooth within a region of similar in size as window w.
 5. Themethod of claim 4, wherein parameter ε determines the smoothness of thedecomposition into said two color modes as represented by the smoothnessof α estimates.
 6. The method of claim 4, wherein α_(i) is defined as:${\alpha_{i} = {{\sum\limits_{c}\; {a^{c}I_{i}^{c}}} + b}},{\forall{i \in w}}$where α are weights for each color channel and is constant for thewindow w, b is a model parameter, and c is an index over the colorchannels.
 7. The method of claim 6, wherein the α values are determinedaccording to quadratic cost function J(α)=α^(T)Lα.
 8. The method ofclaim 7, wherein laplacian matrix L is an n×n matrix, where n is thetotal number of pixels in input image I whose (ij)^(th) element is givenby$\sum\limits_{k{{({i,j})} \in w_{k}}}\; \left( {\delta_{ij} - {\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}} \right)$where δ_(ij) is the Kronecker delta, μ_(k) is a 3×1 mean vector ofcolors inside the k^(th) k window with both i and j as members, σ_(k) isthe 3×3 covariance matrix, |w_(k)| is the cardinality of the window andI₃ is the identity matrix of size 3×3.
 9. The method of claim 8, whereinthe laplacian matrix L is further decomposed into a diagonal matrix Dand a weight matrix W with the formulation L=D−W, wherein: diagonalmatrix D is defined with terms D_(ii)=#[k|iεw_(k)] at its diagonal andrepresents the cardinality of the number of windows of which the pixel iis a member; individual terms of weight matrix W are given by$W_{ij} = {\sum\limits_{k{{({i,j})} \in w_{k}}}{\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}$where δ_(ij) is the Kronecker delta, μ_(k) is a 3×1 mean vector ofcolors inside the k^(th) k window with both i and j as members, σ_(k) isthe 3×3 covariance matrix, |w_(k)| is the cardinality of the window andI₃ is the identity matrix of size 3×3; and color affinity W_(ij) for twopixels with the same color is a positive quantity varying with thehomogeneity of the local windows containing the pixels i and j, andcolor affinity for pixels with different color is zero.
 10. The methodof claim 9, wherein at the local minima, the solution α* satisfies afirst order optimality condition L^(T)α*=0, and the optimal condition isdefined as L^(T)α*=(D−W)^(T)α*.
 11. The method of claim 9, wherein termsof W_(i,j) are evaluated by counting the contribution of only the centerpixel's local window defined as the window centered about it.
 12. Themethod of claim 11, wherein the local window w is a 3×3 pixel windowrequiring O(w²) computations.
 13. The method of claim 1, furtherincluding: (b) interpolating said input image I to a resolution higherthan its native resolution; (c) combining the results of step (a) withthe interpolated, higher resolution image of step (b).
 14. The method ofclaim 13, wherein step (a) generates detail line information of saidinput image I, and step (c) adds this line information to theinterpolated, higher resolution image.
 15. The method of claim 1,wherein the laplacian matrix L is a zero-sum filter kernel.
 16. A methodof enlarging an input image I, comprising: providing said digital inputimage, I, in physical memory, said input image having color informationin a native color space for each image pixel; providing a processingdevice coupled to said physical memory to access said input image I andcreate an output image, h, by implementing the following steps: (a)applying the following sequence to said input image I,${h_{\sigma,ɛ}(x)} = \frac{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}$where: (i) x and y are pixel locations over an image grid; (ii) Ω_(x) isthe neighborhood induced around the central pixel x; (iii) f_(s)(x,y)measures the spatial affinity between pixels at x and y; (iv) parameterσ controls an amount of spatial blurring; (v) L_(xy) ^(ε) a laplacianmatrix that provides positive color affinity value for two examinedpixels, x and y, that have the same color and provides a zero coloraffinity value for pixels with different color; (vi) parameter ε is aregularization term whose relative weight determines the smoothness ofoutput image h; (vii) laplacian matrix L is defined as a diagonalmatrix, D, and a weight matrix, W, with the formulation L=D−W, wherediagonal matrix D is further defined as D_(ii)=#[k|iεw_(k)] at itsdiagonal, which represents the cardinality of the number of windows ofwhich the pixel i is a member, individual terms of the weight matrix Ware given by$W_{ij} = {\sum\limits_{k{{({i,j})} \in w_{k}}}{\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}$and weight matrix W_(ij) is evaluated over all possible overlappingwindows that contain the center pixel; and (b) interpolating said inputimage I to a resolution higher than its native resolution; (c) combiningthe results of step (a) with the interpolated, higher resolution imageof step (b).
 17. The method of claim 16, wherein step (b) generatesdetail line information of said input image I, and step (c) adds thisline information to the interpolated, higher resolution image.
 18. Themethod of claim 16, wherein step (a) is applied to said input image I inits native color space.
 19. The method of claim 17, wherein said nativecolor space is the RGB color space.
 20. A method of smoothing an inputimage I, comprising: providing said digital input image, I, in physicalmemory, said input image having color information in a native colorspace for each image pixel; providing a processing device coupled tosaid physical memory to access said input image I and create an outputimage, h, by implementing the following steps: (a) applying thefollowing sequence to said input image I,${h_{\sigma,ɛ}(x)} = \frac{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}{I(y)}} \right)}{\sum\limits_{y \in \Omega_{x}}\; \left( {{f_{s}^{\sigma}\left( {x,y} \right)}L_{xy}^{ɛ}} \right)}$where: (i) x and y are pixel locations over an image grid; (ii) Ω_(x) isthe neighborhood induced around the central pixel x; (iii) f_(s)(x,y)measures the spatial affinity between pixels at x and y; (iv) parameterσ controls an amount of spatial blurring; (v) L_(xy) ^(ε) a laplacianmatrix that provides positive color affinity value for two examinedpixels, x and y, that have the same color and provides a zero coloraffinity value for pixels with different color; (vi) parameter ε is aregularization term whose relative weight determines the smoothness ofoutput image h; (vii) laplacian matrix L is defined as a diagonalmatrix, D, and a weight matrix, W, with the formulation L=D−W, wherediagonal matrix D is further defined as D_(ii)=#[k|iεw_(k)] at itsdiagonal, which represents the cardinality of the number of windows ofwhich the pixel i is a member, individual terms of the weight matrix Ware given by$W_{ij} = {\sum\limits_{k{{({i,j})} \in w_{k}}}{\frac{1}{w_{k}}\left( {1 + {\left( {I_{i} - \mu_{k}} \right)\left( {\sigma_{k} + {\frac{ɛ}{w_{k}}I_{3}}} \right)^{- 1}\left( {I_{j} - \mu_{k}} \right)}} \right)}}$and terms of W_(i,j) are evaluated by counting the contribution of onlythe center pixel's local window defined as the window centered about it.21. The method of claim 20, wherein the local window w is a 3×3 pixelwindow.
 22. The method of claim 20, wherein regularization factor of εis at last 0.1.
 23. The method of claim 20, wherein step (a) is appliedto said input image I in its native color space.
 24. The method of claim23, wherein said native color space is the RGB color space.